##########################################
# Application
# US election
# Substantive interpretation
##########################################
# load library
library(rstan)
# load data
load("replication_AN/cbq_vote_q1.RData")
dat1 <- summary(cbq_fit)
dat1 <- dat1$summary
load("replication_AN/cbq_vote_q2.RData")
dat2 <- summary(cbq_fit)
dat2 <- dat2$summary
load("replication_AN/cbq_vote_q3.RData")
dat3 <- summary(cbq_fit)
dat3 <- dat3$summary
load("replication_AN/cbq_vote_q4.RData")
dat4 <- summary(cbq_fit)
dat4 <- dat4$summary
load("replication_AN/cbq_vote_q5.RData")
dat5 <- summary(cbq_fit)
dat5 <- dat5$summary
load("replication_AN/cbq_vote_q6.RData")
dat6 <- summary(cbq_fit)
dat6 <- dat6$summary
load("replication_AN/cbq_vote_q7.RData")
dat7 <- summary(cbq_fit)
dat7 <- dat7$summary
load("replication_AN/cbq_vote_q8.RData")
dat8 <- summary(cbq_fit)
dat8 <- dat8$summary
load("replication_AN/cbq_vote_q9.RData")
dat9 <- summary(cbq_fit)
dat9 <- dat9$summary
# preprocess
q_mean <- cbind(dat1[,1],
               dat2[,1],
               dat3[,1],
               dat4[,1],
               dat5[,1],
               dat6[,1],
               dat7[,1],
               dat8[,1],
               dat9[,1])
q_mean <- q_mean[-40,]
q_25 <- cbind(dat1[,4],
             dat2[,4],
             dat3[,4],
             dat4[,4],
             dat5[,4],
             dat6[,4],
             dat7[,4],
             dat8[,4],
             dat9[,4])
q_25 <- q_25[-40,]
q_975 <- cbind(dat1[,8],
             dat2[,8],
             dat3[,8],
             dat4[,8],
             dat5[,8],
             dat6[,8],
             dat7[,8],
             dat8[,8],
             dat9[,8])
q_975 <- q_975[-40,]

var_names <- c(
  "Ideological Distance",
  # Bush
  "Felt personal finance were worse",
  "Felt national economy was worse",
  "Oppose government jobs",
  "Oppose government health care",
  "Oppose government minority assistance",
  "Abortion",
  "Region (East)",
  "Region (South)",
  "Region (West)",
  "New or returning voter",
  "Term limits",
  "Felt deficit was a major problem",
  "Democrat",
  "Republican",
  "Gender (Female)",
  "Respondents education",
  "Age: 18-29",
  "Age: 30-44",
  "Age: 45-59",
  # Clinton
  "Felt personal finance were worse",
  "Felt national economy was worse",
  "Oppose government jobs",
  "Oppose government health care",
  "Oppose government minority assistance",
  "Abortion",
  "Region (East)",
  "Region (South)",
  "Region (West)",
  "New or returning voter",
  "Term limits",
  "Felt deficit was a major problem",
  "Democrat",
  "Republican",
  "Gender (Female)",
  "Respondents education",
  "Age: 18-29",
  "Age: 30-44",
  "Age: 45-59"
  ) 

q_sd <- apply(q_mean,1,sd)
q_sd <- cbind(var_names,q_sd)
q_sd1 <- q_sd[-1,]
q_sd2 <- as.data.frame(cbind(q_sd1[1:19,1:2],q_sd1[20:38,2]))
names(q_sd2) <- c("Covariates","Bush","Clinton")
q_sd2$Bush <- as.numeric(as.character(q_sd2$Bush))
q_sd2$Clinton <- as.numeric(as.character(q_sd2$Clinton))

##########################################
# standard deviation of CBQ estimates
function_c <- function(x){
  sprintf("%.3f", x)
}
table2 <- apply(q_sd2[,2:3],1:2,function_c)
q_sd2[,2:3] <- table2
row.names(q_sd2) <- NULL
write.csv(q_sd2, 'tables/table2.csv')

##########################################
# coef plots
var_names0 <- c(
  "Ideological Distance",
  # Bush
  paste(
    "Bush:",
    c("Felt personal finance were worse",
      "Felt national economy was worse",
      "Oppose government jobs",
      "Oppose government health care",
      "Oppose government minority assistance",
      "Abortion",
      "Region (East)",
      "Region (South)",
      "Region (West)",
      "New or returning voter",
      "Term limits",
      "Felt deficit was a major problem",
      "Democrat",
      "Republican",
      "Gender (Female)",
      "Respondents education",
      "Age: 18-29",
      "Age: 30-44",
      "Age: 45-59")
  )
  ,
  # Clinton
  paste(
    "Clinton:",
    c("Felt personal finance were worse",
      "Felt national economy was worse",
      "Oppose government jobs",
      "Oppose government health care",
      "Oppose government minority assistance",
      "Abortion",
      "Region (East)",
      "Region (South)",
      "Region (West)",
      "New or returning voter",
      "Term limits",
      "Felt deficit was a major problem",
      "Democrat",
      "Republican",
      "Gender (Female)",
      "Respondents education",
      "Age: 18-29",
      "Age: 30-44",
      "Age: 45-59")
  )
) 

var_names1 = c(
  "Ideological Distance",
  "Felt personal finance were worse",
  "Felt national economy was worse",
  "Oppose government jobs",
  "Oppose government health care",
  "Oppose government minority assistance",
  "Abortion",
  "Region (East)",
  "Region (South)",
  "Region (West)",
  "New or returning voter",
  "Term limits",
  "Felt deficit was a major problem",
  "Democrat",
  "Republican",
  "Gender (Female)",
  "Respondents education",
  "Age: 18-29",
  "Age: 30-44",
  "Age: 45-59"
) 




###############################################################
# coef plot: Bush

pdf("figures/figure4_appendix.pdf",width=11,height = 14)
par(mfrow=c(5,4),
    cex.lab=1.5)
plot(NA,xlim=c(0,10),ylim=c(min(q_25[1,]),max(q_975[1,])),
     axes=F,
     xlab="Quantiles",
     ylab=""
     ,main = var_names1[1])
axis(1,at=1:9,labels=paste0("0.",1:9))
axis(2)
lines(x=1:9,y=q_mean[1,])
lines(x=1:9,y=q_25[1,],lty="dashed")
lines(x=1:9,y=q_975[1,],lty="dashed")
abline(h=0,lty="dotted",lwd=0.5,col="gray")

plot(NA,xlim=c(0,10),ylim=c(min(q_25[2,]),max(q_975[2,])),
     axes=F,
     xlab="Quantiles",
     ylab=""
     ,main = var_names1[2])
axis(1,at=1:9,labels=paste0("0.",1:9))
axis(2)
lines(x=1:9,y=q_mean[2,])
lines(x=1:9,y=q_25[2,],lty="dashed")
lines(x=1:9,y=q_975[2,],lty="dashed")
abline(h=0,lty="dotted",lwd=0.5,col="gray")

plot(NA,xlim=c(0,10),ylim=c(min(q_25[3,]),max(q_975[3,])),
     axes=F,
     xlab="Quantiles",
     ylab=""
     ,main = var_names1[3])
axis(1,at=1:9,labels=paste0("0.",1:9))
axis(2)
lines(x=1:9,y=q_mean[3,])
lines(x=1:9,y=q_25[3,],lty="dashed")
lines(x=1:9,y=q_975[3,],lty="dashed")
abline(h=0,lty="dotted",lwd=0.5,col="gray")

plot(NA,xlim=c(0,10),ylim=c(min(q_25[4,]),max(q_975[4,])),
     axes=F,
     xlab="Quantiles",
     ylab=""
     ,main = var_names1[4])
axis(1,at=1:9,labels=paste0("0.",1:9))
axis(2)
lines(x=1:9,y=q_mean[4,])
lines(x=1:9,y=q_25[4,],lty="dashed")
lines(x=1:9,y=q_975[4,],lty="dashed")
abline(h=0,lty="dotted",lwd=0.5,col="gray")

plot(NA,xlim=c(0,10),ylim=c(min(q_25[5,]),max(q_975[5,])),
     axes=F,
     xlab="Quantiles",
     ylab=""
     ,main = var_names1[5])
axis(1,at=1:9,labels=paste0("0.",1:9))
axis(2)
lines(x=1:9,y=q_mean[5,])
lines(x=1:9,y=q_25[5,],lty="dashed")
lines(x=1:9,y=q_975[5,],lty="dashed")
abline(h=0,lty="dotted",lwd=0.5,col="gray")

plot(NA,xlim=c(0,10),ylim=c(min(q_25[6,]),max(q_975[6,])),
     axes=F,
     xlab="Quantiles",
     ylab=""
     ,main = var_names1[6])
axis(1,at=1:9,labels=paste0("0.",1:9))
axis(2)
lines(x=1:9,y=q_mean[6,])
lines(x=1:9,y=q_25[6,],lty="dashed")
lines(x=1:9,y=q_975[6,],lty="dashed")
abline(h=0,lty="dotted",lwd=0.5,col="gray")

plot(NA,xlim=c(0,10),ylim=c(min(q_25[7,]),max(q_975[7,])),
     axes=F,
     xlab="Quantiles",
     ylab=""
     ,main = var_names1[7])
axis(1,at=1:9,labels=paste0("0.",1:9))
axis(2)
lines(x=1:9,y=q_mean[7,])
lines(x=1:9,y=q_25[7,],lty="dashed")
lines(x=1:9,y=q_975[7,],lty="dashed")
abline(h=0,lty="dotted",lwd=0.5,col="gray")

plot(NA,xlim=c(0,10),ylim=c(min(q_25[8,]),max(q_975[8,])),
     axes=F,
     xlab="Quantiles",
     ylab=""
     ,main = var_names1[8])
axis(1,at=1:9,labels=paste0("0.",1:9))
axis(2)
lines(x=1:9,y=q_mean[8,])
lines(x=1:9,y=q_25[8,],lty="dashed")
lines(x=1:9,y=q_975[8,],lty="dashed")
abline(h=0,lty="dotted",lwd=0.5,col="gray")

plot(NA,xlim=c(0,10),ylim=c(min(q_25[9,]),max(q_975[9,])),
     axes=F,
     xlab="Quantiles",
     ylab=""
     ,main = var_names1[9])
axis(1,at=1:9,labels=paste0("0.",1:9))
axis(2)
lines(x=1:9,y=q_mean[9,])
lines(x=1:9,y=q_25[9,],lty="dashed")
lines(x=1:9,y=q_975[9,],lty="dashed")
abline(h=0,lty="dotted",lwd=0.5,col="gray")

plot(NA,xlim=c(0,10),ylim=c(min(q_25[10,]),max(q_975[10,])),
     axes=F,
     xlab="Quantiles",
     ylab=""
     ,main = var_names1[10])
axis(1,at=1:9,labels=paste0("0.",1:9))
axis(2)
lines(x=1:9,y=q_mean[10,])
lines(x=1:9,y=q_25[10,],lty="dashed")
lines(x=1:9,y=q_975[10,],lty="dashed")
abline(h=0,lty="dotted",lwd=0.5,col="gray")

plot(NA,xlim=c(0,10),ylim=c(min(q_25[11,]),max(q_975[11,])),
     axes=F,
     xlab="Quantiles",
     ylab=""
     ,main = var_names1[11])
axis(1,at=1:9,labels=paste0("0.",1:9))
axis(2)
lines(x=1:9,y=q_mean[11,])
lines(x=1:9,y=q_25[11,],lty="dashed")
lines(x=1:9,y=q_975[11,],lty="dashed")
abline(h=0,lty="dotted",lwd=0.5,col="gray")

plot(NA,xlim=c(0,10),ylim=c(min(q_25[12,]),max(q_975[12,])),
     axes=F,
     xlab="Quantiles",
     ylab=""
     ,main = var_names1[12])
axis(1,at=1:9,labels=paste0("0.",1:9))
axis(2)
lines(x=1:9,y=q_mean[12,])
lines(x=1:9,y=q_25[12,],lty="dashed")
lines(x=1:9,y=q_975[12,],lty="dashed")
abline(h=0,lty="dotted",lwd=0.5,col="gray")

plot(NA,xlim=c(0,10),ylim=c(min(q_25[13,]),max(q_975[13,])),
     axes=F,
     xlab="Quantiles",
     ylab=""
     ,main = var_names1[13])
axis(1,at=1:9,labels=paste0("0.",1:9))
axis(2)
lines(x=1:9,y=q_mean[13,])
lines(x=1:9,y=q_25[13,],lty="dashed")
lines(x=1:9,y=q_975[13,],lty="dashed")
abline(h=0,lty="dotted",lwd=0.5,col="gray")

plot(NA,xlim=c(0,10),ylim=c(min(q_25[14,]),max(q_975[14,])),
     axes=F,
     xlab="Quantiles",
     ylab=""
     ,main = var_names1[14])
axis(1,at=1:9,labels=paste0("0.",1:9))
axis(2)
lines(x=1:9,y=q_mean[14,])
lines(x=1:9,y=q_25[14,],lty="dashed")
lines(x=1:9,y=q_975[14,],lty="dashed")
abline(h=0,lty="dotted",lwd=0.5,col="gray")


plot(NA,xlim=c(0,10),ylim=c(min(q_25[15,]),max(q_975[15,])),
     axes=F,
     xlab="Quantiles",
     ylab=""
     ,main = var_names1[15])
axis(1,at=1:9,labels=paste0("0.",1:9))
axis(2)
lines(x=1:9,y=q_mean[15,])
lines(x=1:9,y=q_25[15,],lty="dashed")
lines(x=1:9,y=q_975[15,],lty="dashed")
abline(h=0,lty="dotted",lwd=0.5,col="gray")

plot(NA,xlim=c(0,10),ylim=c(min(q_25[16,]),max(q_975[16,])),
     axes=F,
     xlab="Quantiles",
     ylab=""
     ,main = var_names1[16])
axis(1,at=1:9,labels=paste0("0.",1:9))
axis(2)
lines(x=1:9,y=q_mean[16,])
lines(x=1:9,y=q_25[16,],lty="dashed")
lines(x=1:9,y=q_975[16,],lty="dashed")
abline(h=0,lty="dotted",lwd=0.5,col="gray")

plot(NA,xlim=c(0,10),ylim=c(min(q_25[17,]),max(q_975[17,])),
     axes=F,
     xlab="Quantiles",
     ylab=""
     ,main = var_names1[17])
axis(1,at=1:9,labels=paste0("0.",1:9))
axis(2)
lines(x=1:9,y=q_mean[17,])
lines(x=1:9,y=q_25[17,],lty="dashed")
lines(x=1:9,y=q_975[17,],lty="dashed")
abline(h=0,lty="dotted",lwd=0.5,col="gray")

plot(NA,xlim=c(0,10),ylim=c(min(q_25[18,]),max(q_975[18,])),
     axes=F,
     xlab="Quantiles",
     ylab=""
     ,main = var_names1[18])
axis(1,at=1:9,labels=paste0("0.",1:9))
axis(2)
lines(x=1:9,y=q_mean[18,])
lines(x=1:9,y=q_25[18,],lty="dashed")
lines(x=1:9,y=q_975[18,],lty="dashed")
abline(h=0,lty="dotted",lwd=0.5,col="gray")

plot(NA,xlim=c(0,10),ylim=c(min(q_25[19,]),max(q_975[19,])),
     axes=F,
     xlab="Quantiles",
     ylab=""
     ,main = var_names1[19])
axis(1,at=1:9,labels=paste0("0.",1:9))
axis(2)
lines(x=1:9,y=q_mean[19,])
lines(x=1:9,y=q_25[19,],lty="dashed")
lines(x=1:9,y=q_975[19,],lty="dashed")
abline(h=0,lty="dotted",lwd=0.5,col="gray")

plot(NA,xlim=c(0,10),ylim=c(min(q_25[20,]),max(q_975[20,])),
     axes=F,
     xlab="Quantiles",
     ylab=""
     ,main = var_names1[20])
axis(1,at=1:9,labels=paste0("0.",1:9))
axis(2)
lines(x=1:9,y=q_mean[20,])
lines(x=1:9,y=q_25[20,],lty="dashed")
lines(x=1:9,y=q_975[20,],lty="dashed")
abline(h=0,lty="dotted",lwd=0.5,col="gray")

dev.off()




##########################################################
# coef plot: Clinton

pdf("figures/figure5_appendix.pdf",width=11,height = 14)
par(mfrow=c(5,4),
    cex.lab=1.5)
plot(NA,xlim=c(0,10),ylim=c(min(q_25[1,]),max(q_975[1,])),
     axes=F,
     xlab="Quantiles",
     ylab=""
     ,main = var_names1[1])
axis(1,at=1:9,labels=paste0("0.",1:9))
axis(2)
lines(x=1:9,y=q_mean[1,])
lines(x=1:9,y=q_25[1,],lty="dashed")
lines(x=1:9,y=q_975[1,],lty="dashed")
abline(h=0,lty="dotted",lwd=0.5,col="gray")

plot(NA,xlim=c(0,10),ylim=c(min(q_25[21,]),max(q_975[21,])),
     axes=F,
     xlab="Quantiles",
     ylab=""
     ,main = var_names1[2])
axis(1,at=1:9,labels=paste0("0.",1:9))
axis(2)
lines(x=1:9,y=q_mean[21,])
lines(x=1:9,y=q_25[21,],lty="dashed")
lines(x=1:9,y=q_975[21,],lty="dashed")
abline(h=0,lty="dotted",lwd=0.5,col="gray")


plot(NA,xlim=c(0,10),ylim=c(min(q_25[22,]),max(q_975[22,])),
     axes=F,
     xlab="Quantiles",
     ylab=""
     ,main = var_names1[3])
axis(1,at=1:9,labels=paste0("0.",1:9))
axis(2)
lines(x=1:9,y=q_mean[22,])
lines(x=1:9,y=q_25[22,],lty="dashed")
lines(x=1:9,y=q_975[22,],lty="dashed")
abline(h=0,lty="dotted",lwd=0.5,col="gray")

plot(NA,xlim=c(0,10),ylim=c(min(q_25[23,]),max(q_975[23,])),
     axes=F,
     xlab="Quantiles",
     ylab=""
     ,main = var_names1[4])
axis(1,at=1:9,labels=paste0("0.",1:9))
axis(2)
lines(x=1:9,y=q_mean[23,])
lines(x=1:9,y=q_25[23,],lty="dashed")
lines(x=1:9,y=q_975[23,],lty="dashed")
abline(h=0,lty="dotted",lwd=0.5,col="gray")

plot(NA,xlim=c(0,10),ylim=c(min(q_25[24,]),max(q_975[24,])),
     axes=F,
     xlab="Quantiles",
     ylab=""
     ,main = var_names1[5])
axis(1,at=1:9,labels=paste0("0.",1:9))
axis(2)
lines(x=1:9,y=q_mean[24,])
lines(x=1:9,y=q_25[24,],lty="dashed")
lines(x=1:9,y=q_975[24,],lty="dashed")
abline(h=0,lty="dotted",lwd=0.5,col="gray")

plot(NA,xlim=c(0,10),ylim=c(min(q_25[25,]),max(q_975[25,])),
     axes=F,
     xlab="Quantiles",
     ylab=""
     ,main = var_names1[6])
axis(1,at=1:9,labels=paste0("0.",1:9))
axis(2)
lines(x=1:9,y=q_mean[25,])
lines(x=1:9,y=q_25[25,],lty="dashed")
lines(x=1:9,y=q_975[25,],lty="dashed")
abline(h=0,lty="dotted",lwd=0.5,col="gray")

plot(NA,xlim=c(0,10),ylim=c(min(q_25[26,]),max(q_975[26,])),
     axes=F,
     xlab="Quantiles",
     ylab=""
     ,main = var_names1[7])
axis(1,at=1:9,labels=paste0("0.",1:9))
axis(2)
lines(x=1:9,y=q_mean[26,])
lines(x=1:9,y=q_25[26,],lty="dashed")
lines(x=1:9,y=q_975[26,],lty="dashed")
abline(h=0,lty="dotted",lwd=0.5,col="gray")

plot(NA,xlim=c(0,10),ylim=c(min(q_25[27,]),max(q_975[27,])),
     axes=F,
     xlab="Quantiles",
     ylab=""
     ,main = var_names1[8])
axis(1,at=1:9,labels=paste0("0.",1:9))
axis(2)
lines(x=1:9,y=q_mean[27,])
lines(x=1:9,y=q_25[27,],lty="dashed")
lines(x=1:9,y=q_975[27,],lty="dashed")
abline(h=0,lty="dotted",lwd=0.5,col="gray")

plot(NA,xlim=c(0,10),ylim=c(min(q_25[28,]),max(q_975[28,])),
     axes=F,
     xlab="Quantiles",
     ylab=""
     ,main = var_names1[9])
axis(1,at=1:9,labels=paste0("0.",1:9))
axis(2)
lines(x=1:9,y=q_mean[28,])
lines(x=1:9,y=q_25[28,],lty="dashed")
lines(x=1:9,y=q_975[28,],lty="dashed")
abline(h=0,lty="dotted",lwd=0.5,col="gray")

plot(NA,xlim=c(0,10),ylim=c(min(q_25[29,]),max(q_975[29,])),
     axes=F,
     xlab="Quantiles",
     ylab=""
     ,main = var_names1[10])
axis(1,at=1:9,labels=paste0("0.",1:9))
axis(2)
lines(x=1:9,y=q_mean[29,])
lines(x=1:9,y=q_25[29,],lty="dashed")
lines(x=1:9,y=q_975[29,],lty="dashed")
abline(h=0,lty="dotted",lwd=0.5,col="gray")


plot(NA,xlim=c(0,10),ylim=c(min(q_25[30,]),max(q_975[30,])),
     axes=F,
     xlab="Quantiles",
     ylab=""
     ,main = var_names1[11])
axis(1,at=1:9,labels=paste0("0.",1:9))
axis(2)
lines(x=1:9,y=q_mean[30,])
lines(x=1:9,y=q_25[30,],lty="dashed")
lines(x=1:9,y=q_975[30,],lty="dashed")
abline(h=0,lty="dotted",lwd=0.5,col="gray")

plot(NA,xlim=c(0,10),ylim=c(min(q_25[31,]),max(q_975[31,])),
     axes=F,
     xlab="Quantiles",
     ylab=""
     ,main = var_names1[12])
axis(1,at=1:9,labels=paste0("0.",1:9))
axis(2)
lines(x=1:9,y=q_mean[31,])
lines(x=1:9,y=q_25[31,],lty="dashed")
lines(x=1:9,y=q_975[31,],lty="dashed")
abline(h=0,lty="dotted",lwd=0.5,col="gray")

plot(NA,xlim=c(0,10),ylim=c(min(q_25[32,]),max(q_975[32,])),
     axes=F,
     xlab="Quantiles",
     ylab=""
     ,main = var_names1[13])
axis(1,at=1:9,labels=paste0("0.",1:9))
axis(2)
lines(x=1:9,y=q_mean[32,])
lines(x=1:9,y=q_25[32,],lty="dashed")
lines(x=1:9,y=q_975[32,],lty="dashed")
abline(h=0,lty="dotted",lwd=0.5,col="gray")

plot(NA,xlim=c(0,10),ylim=c(min(q_25[33,]),max(q_975[33,])),
     axes=F,
     xlab="Quantiles",
     ylab=""
     ,main = var_names1[14])
axis(1,at=1:9,labels=paste0("0.",1:9))
axis(2)
lines(x=1:9,y=q_mean[33,])
lines(x=1:9,y=q_25[33,],lty="dashed")
lines(x=1:9,y=q_975[33,],lty="dashed")
abline(h=0,lty="dotted",lwd=0.5,col="gray")

plot(NA,xlim=c(0,10),ylim=c(min(q_25[34,]),max(q_975[34,])),
     axes=F,
     xlab="Quantiles",
     ylab=""
     ,main = var_names1[15])
axis(1,at=1:9,labels=paste0("0.",1:9))
axis(2)
lines(x=1:9,y=q_mean[34,])
lines(x=1:9,y=q_25[34,],lty="dashed")
lines(x=1:9,y=q_975[34,],lty="dashed")
abline(h=0,lty="dotted",lwd=0.5,col="gray")

plot(NA,xlim=c(0,10),ylim=c(min(q_25[35,]),max(q_975[35,])),
     axes=F,
     xlab="Quantiles",
     ylab=""
     ,main = var_names1[16])
axis(1,at=1:9,labels=paste0("0.",1:9))
axis(2)
lines(x=1:9,y=q_mean[35,])
lines(x=1:9,y=q_25[35,],lty="dashed")
lines(x=1:9,y=q_975[35,],lty="dashed")
abline(h=0,lty="dotted",lwd=0.5,col="gray")

plot(NA,xlim=c(0,10),ylim=c(min(q_25[36,]),max(q_975[36,])),
     axes=F,
     xlab="Quantiles",
     ylab=""
     ,main = var_names1[17])
axis(1,at=1:9,labels=paste0("0.",1:9))
axis(2)
lines(x=1:9,y=q_mean[36,])
lines(x=1:9,y=q_25[36,],lty="dashed")
lines(x=1:9,y=q_975[36,],lty="dashed")
abline(h=0,lty="dotted",lwd=0.5,col="gray")

plot(NA,xlim=c(0,10),ylim=c(min(q_25[37,]),max(q_975[37,])),
     axes=F,
     xlab="Quantiles",
     ylab=""
     ,main = var_names1[18])
axis(1,at=1:9,labels=paste0("0.",1:9))
axis(2)
lines(x=1:9,y=q_mean[37,])
lines(x=1:9,y=q_25[37,],lty="dashed")
lines(x=1:9,y=q_975[37,],lty="dashed")
abline(h=0,lty="dotted",lwd=0.5,col="gray")

plot(NA,xlim=c(0,10),ylim=c(min(q_25[38,]),max(q_975[38,])),
     axes=F,
     xlab="Quantiles",
     ylab=""
     ,main = var_names1[19])
axis(1,at=1:9,labels=paste0("0.",1:9))
axis(2)
lines(x=1:9,y=q_mean[38,])
lines(x=1:9,y=q_25[38,],lty="dashed")
lines(x=1:9,y=q_975[38,],lty="dashed")
abline(h=0,lty="dotted",lwd=0.5,col="gray")

plot(NA,xlim=c(0,10),ylim=c(min(q_25[39,]),max(q_975[39,])),
     axes=F,
     xlab="Quantiles",
     ylab=""
     ,main = var_names1[20])
axis(1,at=1:9,labels=paste0("0.",1:9))
axis(2)
lines(x=1:9,y=q_mean[39,])
lines(x=1:9,y=q_25[39,],lty="dashed")
lines(x=1:9,y=q_975[39,],lty="dashed")
abline(h=0,lty="dotted",lwd=0.5,col="gray")


dev.off()


##################################
# read in dataset

load("replication_AN/vote_dat_origin.RData")
load("replication_AN/vote_dat_sim.RData")
source("replication_AN/cbq_function.R")


#############################################################
# calculate predicted probabilities
x0 <- as.matrix(subset(vote_dat,select = c(
  dist , respfinp.bush , natlec.bush , respgjob.bush , resphlth.bush , 
  respblk.bush , respab.bush , east.bush , south.bush , west.bush , newvoter.bush , termlim.bush ,
  deficit.bush , dem.bush , rep.bush , women.bush , educ.bush , age1829.bush , age3044.bush , age4559.bush
  , respfinp.clinton , natlec.clinton , respgjob.clinton , resphlth.clinton , 
  respblk.clinton , respab.clinton , east.clinton , south.clinton , west.clinton , newvoter.clinton , termlim.clinton ,
  deficit.clinton , dem.clinton , rep.clinton , women.clinton , educ.clinton , age1829.clinton , age3044.clinton , age4559.clinton
)))
  
xb0 <- x0 %*% q_mean

pred_prob1 <- pald(xb0[,1],0.1)
pred_prob2 <- pald(xb0[,2],0.2)
pred_prob3 <- pald(xb0[,3],0.3)
pred_prob4 <- pald(xb0[,4],0.4)
pred_prob5 <- pald(xb0[,5],0.5)
pred_prob6 <- pald(xb0[,6],0.6)
pred_prob7 <- pald(xb0[,7],0.7)
pred_prob8 <- pald(xb0[,8],0.8)
pred_prob9 <- pald(xb0[,9],0.9)
tmp <- cbind(vote_dat$chid,pred_prob1,pred_prob2,pred_prob3,pred_prob4,pred_prob5,pred_prob6,pred_prob7,
            pred_prob8,pred_prob9,vote_dat$preschc)

pred1 <- rep(0,909)
for (i in 1:909){
  nnn <- tmp[which(tmp[,1]==i),2]
  pred1[i] <- which.max(nnn)
}

pred2 <- rep(0,909)
for (i in 1:909){
  nnn <- tmp[which(tmp[,1]==i),3]
  pred2[i] <- which.max(nnn)
}

pred3 <- rep(0,909)
for (i in 1:909){
  nnn <- tmp[which(tmp[,1]==i),4]
  pred3[i] <- which.max(nnn)
}

pred4 <- rep(0,909)
for (i in 1:909){
  nnn <- tmp[which(tmp[,1]==i),5]
  pred4[i] <- which.max(nnn)
}

pred5 <- rep(0,909)
for (i in 1:909){
  nnn <- tmp[which(tmp[,1]==i),6]
  pred5[i] <- which.max(nnn)
}

pred6 <- rep(0,909)
for (i in 1:909){
  nnn <- tmp[which(tmp[,1]==i),7]
  pred6[i] <- which.max(nnn)
}

pred7 <- rep(0,909)
for (i in 1:909){
  nnn <- tmp[which(tmp[,1]==i),8]
  pred7[i] <- which.max(nnn)
}

pred8 <- rep(0,909)
for (i in 1:909){
  nnn <- tmp[which(tmp[,1]==i),9]
  pred8[i] <- which.max(nnn)
}

pred9 <- rep(0,909)
for (i in 1:909){
  nnn <- tmp[which(tmp[,1]==i),10]
  pred9[i] <- which.max(nnn)
}

# load original dataset
dat <- read.table("replication_AN/nes9212r.asc")
varnames <- c("preschc", # choice: 1 Bush 2 Clinton 3 Perot
             "educ", # Respondents education
             "east", # Region (East)
             "south", # Region (South)
             "west", # Region (West)
             "women", # Gender (Female)
             "respfinp", # Felt personal finance were worse
             "natlec", # Felt national economy was worse
             "bclibdis", # Ideological Distance
             "gblibdis", # Ideological Distance
             "rplibdis", # Ideological Distance
             "dem", # Democrat
             "rep", # Republican
             "respgjob", # Oppose government jobs
             "resphlth", # Oppose government health care
             "respblk", # Oppose government minority assistance
             "respab", # Abortion
             "termlim", # Term limits
             "age1829", # Age: 18-29
             "age3044", # Age: 30-44
             "age4559", # Age: 45-59
             "newvoter", # New or returning voter
             "deficit") # Felt deficit was a major problem

names(dat) <- varnames
dat$choice <- "Bush"
dat$choice[which(dat$preschc==2)] <- "Clinton"
dat$choice[which(dat$preschc==3)] <- "Perot"

names(dat)[9] <- "dist.Clinton"
names(dat)[10] <- "dist.Bush"
names(dat)[11] <- "dist.Perot"

real_choice <- dat$preschc

# original prediction accuracy: 70.6%
pred_rate1 <- length(which(pred1==real_choice))/909
pred_rate2 <- length(which(pred2==real_choice))/909
pred_rate3 <- length(which(pred3==real_choice))/909
pred_rate4 <- length(which(pred4==real_choice))/909
pred_rate5 <- length(which(pred5==real_choice))/909
pred_rate6 <- length(which(pred6==real_choice))/909
pred_rate7 <- length(which(pred7==real_choice))/909
pred_rate8 <- length(which(pred8==real_choice))/909
pred_rate9 <- length(which(pred9==real_choice))/909

pred_rates <- sapply(list(pred_rate1,
            pred_rate2,
            pred_rate3,
            pred_rate4,
            pred_rate5,
            pred_rate6,
            pred_rate7,
            pred_rate8,
            pred_rate9),
       round,3)
min(pred_rates) # 74.4%
max(pred_rates) # 75.6%
# compared to mprobit: 74.1%
##########################################################

pos_seq <- seq(from = 0.1, to = 7, by = 0.02)
preds1 <- matrix(1,ncol = 3,nrow = 909)
preds2 <- matrix(1,ncol = 3,nrow = 909)
preds3 <- matrix(1,ncol = 3,nrow = 909)
preds4 <- matrix(1,ncol = 3,nrow = 909)
preds5 <- matrix(1,ncol = 3,nrow = 909)
preds6 <- matrix(1,ncol = 3,nrow = 909)
preds7 <- matrix(1,ncol = 3,nrow = 909)
preds8 <- matrix(1,ncol = 3,nrow = 909)
preds9 <- matrix(1,ncol = 3,nrow = 909)

for (j in 1:length(pos_seq)){
  vote_dat2$dist[which(vote_dat2$bush==1)] <- (vote_dat2$resplib[which(vote_dat2$bush==1)] - pos_seq[j])^2
  x1 <- as.matrix(subset(vote_dat2,select = c(
    dist , respfinp.bush , natlec.bush , respgjob.bush , resphlth.bush , 
    respblk.bush , respab.bush , east.bush , south.bush , west.bush , newvoter.bush , termlim.bush ,
    deficit.bush , dem.bush , rep.bush , women.bush , educ.bush , age1829.bush , age3044.bush , age4559.bush
    , respfinp.clinton , natlec.clinton , respgjob.clinton , resphlth.clinton , 
    respblk.clinton , respab.clinton , east.clinton , south.clinton , west.clinton , newvoter.clinton , termlim.clinton ,
    deficit.clinton , dem.clinton , rep.clinton , women.clinton , educ.clinton , age1829.clinton , age3044.clinton , age4559.clinton
  )))
  xb1 <- x1 %*% q_mean
  
  pred0_prob1 <- pald(xb1[,1],0.1)
  pred0_prob2 <- pald(xb1[,2],0.2)
  pred0_prob3 <- pald(xb1[,3],0.3)
  pred0_prob4 <- pald(xb1[,4],0.4)
  pred0_prob5 <- pald(xb1[,5],0.5)
  pred0_prob6 <- pald(xb1[,6],0.6)
  pred0_prob7 <- pald(xb1[,7],0.7)
  pred0_prob8 <- pald(xb1[,8],0.8)
  pred0_prob9 <- pald(xb1[,9],0.9)
  
  tmp <- cbind(vote_dat2$chid,
              pred0_prob1,
              pred0_prob2,
              pred0_prob3,
              pred0_prob4,
              pred0_prob5,
              pred0_prob6,
              pred0_prob7,
              pred0_prob8,
              pred0_prob9,
              vote_dat2$preschc)
  
  pred1 <- rep(0,909)
  for (i in 1:909){
    nnn <- tmp[which(tmp[,1]==i),2]
    pred1[i] <- which.max(nnn)
  }
  
  pred2 <- rep(0,909)
  for (i in 1:909){
    nnn <- tmp[which(tmp[,1]==i),3]
    pred2[i] <- which.max(nnn)
  }
  
  pred3 <- rep(0,909)
  for (i in 1:909){
    nnn <- tmp[which(tmp[,1]==i),4]
    pred3[i] <- which.max(nnn)
  }
  
  pred4 <- rep(0,909)
  for (i in 1:909){
    nnn <- tmp[which(tmp[,1]==i),5]
    pred4[i] <- which.max(nnn)
  }
  
  pred5 <- rep(0,909)
  for (i in 1:909){
    nnn <- tmp[which(tmp[,1]==i),6]
    pred5[i] <- which.max(nnn)
  }
  
  pred6 <- rep(0,909)
  for (i in 1:909){
    nnn <- tmp[which(tmp[,1]==i),7]
    pred6[i] <- which.max(nnn)
  }
  
  pred7 <- rep(0,909)
  for (i in 1:909){
    nnn <- tmp[which(tmp[,1]==i),8]
    pred7[i] <- which.max(nnn)
  }
  
  pred8 <- rep(0,909)
  for (i in 1:909){
    nnn <- tmp[which(tmp[,1]==i),9]
    pred8[i] <- which.max(nnn)
  }
  
  pred9 <- rep(0,909)
  for (i in 1:909){
    nnn <- tmp[which(tmp[,1]==i),10]
    pred9[i] <- which.max(nnn)
  }
  
  preds1[j,1] <- length(which(pred1==1))/909
  preds1[j,2] <- length(which(pred1==2))/909
  preds1[j,3] <- length(which(pred1==3))/909
  
  preds2[j,1] <- length(which(pred2==1))/909
  preds2[j,2] <- length(which(pred2==2))/909
  preds2[j,3] <- length(which(pred2==3))/909
  
  preds3[j,1] <- length(which(pred3==1))/909
  preds3[j,2] <- length(which(pred3==2))/909
  preds3[j,3] <- length(which(pred3==3))/909
  
  preds4[j,1] <- length(which(pred4==1))/909
  preds4[j,2] <- length(which(pred4==2))/909
  preds4[j,3] <- length(which(pred4==3))/909
  
  preds5[j,1] <- length(which(pred5==1))/909
  preds5[j,2] <- length(which(pred5==2))/909
  preds5[j,3] <- length(which(pred5==3))/909
  
  preds6[j,1] <- length(which(pred6==1))/909
  preds6[j,2] <- length(which(pred6==2))/909
  preds6[j,3] <- length(which(pred6==3))/909
  
  preds7[j,1] <- length(which(pred7==1))/909
  preds7[j,2] <- length(which(pred7==2))/909
  preds7[j,3] <- length(which(pred7==3))/909
  
  preds8[j,1] <- length(which(pred8==1))/909
  preds8[j,2] <- length(which(pred8==2))/909
  preds8[j,3] <- length(which(pred8==3))/909
  
  preds9[j,1] <- length(which(pred9==1))/909
  preds9[j,2] <- length(which(pred9==2))/909
  preds9[j,3] <- length(which(pred9==3))/909
  
}

preds_all_q <- list(preds1,
                   preds2,
                   preds3,
                   preds4,
                   preds5,
                   preds6,
                   preds7,
                   preds8,
                   preds9)
save(preds_all_q,file="replication_AN/pred_vote_share_all_quantile.RData")
# 1: Bush 2: Clinton 3: Perot


# load original dataset for plot
dat <- read.table("replication_AN/nes9212r.asc")
varnames <- c("preschc", # choice: 1 Bush 2 Clinton 3 Perot
             "educ", # Respondents education
             "east", # Region (East)
             "south", # Region (South)
             "west", # Region (West)
             "women", # Gender (Female)
             "respfinp", # Felt personal finance were worse
             "natlec", # Felt national economy was worse
             "bclibdis", # Ideological Distance
             "gblibdis", # Ideological Distance
             "rplibdis", # Ideological Distance
             "dem", # Democrat
             "rep", # Republican
             "respgjob", # Oppose government jobs
             "resphlth", # Oppose government health care
             "respblk", # Oppose government minority assistance
             "respab", # Abortion
             "termlim", # Term limits
             "age1829", # Age: 18-29
             "age3044", # Age: 30-44
             "age4559", # Age: 45-59
             "newvoter", # New or returning voter
             "deficit") # Felt deficit was a major problem

names(dat) <- varnames
dat$choice <- "Bush"
dat$choice[which(dat$preschc==2)] <- "Clinton"
dat$choice[which(dat$preschc==3)] <- "Perot"

names(dat)[9] <- "dist.Clinton"
names(dat)[10] <- "dist.Bush"
names(dat)[11] <- "dist.Perot"

load("replication_AN/pred_vote_share_all_quantile.RData")
pos_seq <- seq(from = 0.1, to = 7, by = 0.02)
preds1 <- preds_all_q[[1]]
preds2 <- preds_all_q[[2]]
preds3 <- preds_all_q[[3]]
preds4 <- preds_all_q[[4]]
preds5 <- preds_all_q[[5]]
preds6 <- preds_all_q[[6]]
preds7 <- preds_all_q[[7]]
preds8 <- preds_all_q[[8]]
preds9 <- preds_all_q[[9]]

pdf("figures/figure4.pdf",width = 10,height = 7)
plot(NA,xlim = c(0,7.1),ylim=c(0,0.6),
     xlab = "Bush's position on 7-point scale",
     ylab="Vote share",
     bty="l")
lines(pos_seq,preds1[1:346,1],lty="dotted") # Perot
lines(pos_seq,preds1[1:346,2],lty="dashed") # Bush
lines(pos_seq,preds1[1:346,3],lty="solid") # Clinton

lines(pos_seq,preds2[1:346,1],lty="dotted") # Perot
lines(pos_seq,preds2[1:346,2],lty="dashed") # Bush
lines(pos_seq,preds2[1:346,3],lty="solid") # Clinton

lines(pos_seq,preds3[1:346,1],lty="dotted") # Perot
lines(pos_seq,preds3[1:346,2],lty="dashed") # Bush
lines(pos_seq,preds3[1:346,3],lty="solid") # Clinton

lines(pos_seq,preds4[1:346,1],lty="dotted") # Perot
lines(pos_seq,preds4[1:346,2],lty="dashed") # Bush
lines(pos_seq,preds4[1:346,3],lty="solid") # Clinton

lines(pos_seq,preds5[1:346,1],lty="dotted") # Perot
lines(pos_seq,preds5[1:346,2],lty="dashed") # Bush
lines(pos_seq,preds5[1:346,3],lty="solid") # Clinton

lines(pos_seq,preds6[1:346,1],lty="dotted") # Perot
lines(pos_seq,preds6[1:346,2],lty="dashed") # Bush
lines(pos_seq,preds6[1:346,3],lty="solid") # Clinton

lines(pos_seq,preds7[1:346,1],lty="dotted") # Perot
lines(pos_seq,preds7[1:346,2],lty="dashed") # Bush
lines(pos_seq,preds7[1:346,3],lty="solid") # Clinton

lines(pos_seq,preds8[1:346,1],lty="dotted") # Perot
lines(pos_seq,preds8[1:346,2],lty="dashed") # Bush
lines(pos_seq,preds8[1:346,3],lty="solid") # Clinton

lines(pos_seq,preds9[1:346,1],lty="dotted") # Perot
lines(pos_seq,preds9[1:346,2],lty="dashed") # Bush
lines(pos_seq,preds9[1:346,3],lty="solid") # Clinton

load("replication_AN/pred_voteshare.RData")

lines(pos_seq,preds[2:347,1],lty="dotted",col="red",lwd=2) # Perot
lines(pos_seq,preds[2:347,2],lty="dashed",col="red",lwd=2) # Bush
lines(pos_seq,preds[2:347,3],lty="solid",col="red",lwd=2) # Clinton

text(rep(5.32,3)-c(0.6),table(dat$preschc)/909,c("Bush","Clinton","Perot"))


points(rep(5.32,3),table(dat$preschc)/909,col="red",pch=15)
legend("bottomright",legend=c("Clinton","Bush","Perot"),lty=c("dashed","dotted","solid"),
       title = "Predicted vote share",
       bty="n",
       ncol = 3)

legend("bottomleft",legend=c("CBQ","Multinomial probit"),lty=c("solid","solid"),
       col=c("black","red"),
       title = "Model",
       bty="n",
       ncol = 2)
dev.off()

